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Abstract: A 2+1 dimensional fermion field theory is proposed as a model for the low- 
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component Dirac fermions moving in the plane and interacting via a contact interaction 
between charge densities. For strong couplings there is a continuous transition to a Mott 
insulting phase. We present results of an extensive numerical study of the model's critical 
region, including the order parameter, its associated susceptibility, and for the first time 
the quasiparticle propagator. The data enables an extraction of the critical exponents 
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from our value of the critical coupling. 
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1 Introduction 

There has been much recent interest in the remarkable electronic properties of graphene 
[Ij (see also a recent, comprehensive review in [2]). It appears that they arise from the 
low-energy spectrum of excitations being equivalent to that of a two-dimensional gas 
of relativistic fermions, with the case of undoped (ie. neutral) graphene corresponding 
to zero net particle number. In brief, for a carbon monolayer with one mobile electron 
per atom, a simple tight-binding model predicts a linear dispersion relation centred on 
zeroes located at the six corners of the first Brillouin zone. It is possible to rewrite 
the Hamiltonian for single-particle excitations in Dirac form with Nf = 2 flavors of 
four-component spinor ip, the counting of degrees of freedom coming from 2 C atoms 
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per unit cell x 2 zeroes per zone x 2 physical spin components per electron. Electron 
propagation within the monolayer is thus relativistic, albeit with speed vp ~ c/300. 

The charge carriers in graphene can be either electrons ("particles") or holes ("an- 
tiparticles") and are characterised by a very high value of mobility /i = a/ne (where 
0" is electrical conductivity and n the carrier density), more than twice that of the 
highest mobility conventional semi-conductor, and several orders of magnitude greater 
than a typical metal at room temperature. This gives graphene the potential to be of 
great technological significance in the construction of fast electronic devices. The naive 
tight-binding model suggests, and experiments with graphene based on a Si02 substrate 
confirm, that graphene remains a conductor (technically a semimetal) for all values of 
the gate voltage, ie. even when the carrier density formally vanishes, because there is 
no gap in the energy spectrum at the Dirac points. The presence of a small gap would, 
however, be extremely valuable for electronics applications, because it would increase 
the effective on-off current flow ratio needed for device stability [3j . 

More sophisticated theoretical approaches to graphene must take inter-electron in- 
teractions into account. In this paper we build upon an approach which treats the 
low energy fermion excitations using a 2+1-dimensional relativistic quantum field the- 
ory [H El |6] . The interaction between electrons is assumed to be an unscreened Coulomb 
potential; because vp ^ c this can be treated as "instantaneous", meaning that the field 
theory is necessarily non-local. The strength of the Coulomb interaction is variable, 
since it depends on the dielectric constant e of the underlying substrate. It can be 
parametrised by an effective fine structure constant 

a = -^^- 0(1), (1) 

so that the problem is strongly-interacting. The possibility then opens up of the dis- 
ruption of the free-field ground state by an excitonic condensate, ie. one formed from 
tightly-bound electron-hole pairs, with the effect of opening up a gap A > at the Dirac 
point and making the ground state of undoped graphene a Mott insulator. 

The analogous phenomenon in particle physics is described using different language: 
we say that the global chiral symmetry of the model, which prevents generation of a 
fermion mass through quantum corrections to all orders in perturbation theory, is spon- 
taneously broken by the formation of a chiral condensate. Dynamical mass generation 
is therefore inherently non-perturbative, and must be addressed either by self-consistent 
analytic methods or, as in this paper, by numerical simulation of a lattice-regularised 
version of the field theory. To date there have been two distinct approaches taken. In a 
series of papers, Drut and Lahde 13 [H [9] have simulated a formulation of the graphene 
field theory based on lattice gauge theory, in which electrostatic degrees of freedom are 
formulated on a 3 + 1-dimensional lattice, while the electron fields are restricted to a 
2+1 dimensional slice. By contrast, our formulation [lU] is entirely 2-1-1 dimensional, 
and is in essence a non-covariant form of the Thirring model. Both numerical calcula- 
tions support the hypothesis proposed in [6] that the semimetal and insulator phases 
are separated by a line of second order phase transitions in the {a, Nf) plane, starting 
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at a point {oo, Nfc) and running in the direction of decreasing a and decreasing Nf. A 
very similar situation pertains in the 2+ld Thirring model [HI [12]. In [10] we studied 
chiral symmetry breaking with variable Nf in the strong-coupling limit and estimated 
Nfc = 4.8(2). Each point on the line with integer Nf < Nfc defines a quantum crit- 
ical point (QCP) whose properties are characterised by a set of iV/-dependent critical 



In the current paper we present an extensive numerical study of the semimetal- 
insulator phase transition for the case Nf = 2 which is of direct physical interest; only 
preliminary results were available in [10]. The graphene model we study will be presented 
in detail, both as a continuum field theory and as a lattice model, in the following section, 
but it is appropriate to preface that with some remarks. Because our approach is based 
upon a local quantum field theory in 2+ld, it is unable to capture the long-range 1/r 
nature of the unscreened Coulomb potential assumed in [U [5l [6]. We have argued in 
[TU] that this is unimportant in the strong-coupling limit a — > cx) where electron-hole 
pair polarisation effects dominate the long-range physics; however for finite a our model 
is in principle different both from the continuum approaches [H El [6] and the lattice 
gauge theory approach of [71IH119]. We do not exclude the possibility, however, that the 
universal behaviour at the QCP remains the same, even for Nf < Nfc- 

There are two main benefits of our approach. Firstly, the simplicity of our model 
and the fact that it is formulated directly on a 2+ld spacetime lattice mean that we 
have been able to perform accurate simulations on a range of system sizes L^x Lt, yield- 
ing control over finite size artifacts and hence access to the model's critical properties. 
Secondly, the fact that our model is not a gauge theory permits a definition of the quasi- 
particle correlation function without any need for gauge fixing, which is known to be a 
major source of statistical noise in similar model systems, eg. [13]. We are able here to 
present the first numerical study of the quasiparticle propagator, which both explicitly 
demonstrates gap generation as the coupling strength g"^ is increased beyond Qc, and 
broadens the scope of the critical analysis; for the first time we are able to present an 
estimate for the dynamical critical exponent z which governs the different scaling of the 
correlation length in spacelike and timelike directions. 

The remainder of the paper is organised as follows. In Sec. [2] we lay out the model to 
be studied in both continuum and lattice formulations, and discuss its relation with other 
models studied in the literature and its applicability to graphene. Our numerical results 
are presented in Sec. [S] Sec. l3.1l focusses on the chiral order parameter and its associated 
susceptibility, and fits data to a renormalisation-group inspired critical equation of state 
yielding estimates for the critical coupling and exponents 6 and /?; Sec. 13.21 presents 
an analysis of the quasiparticle propagator and shows how both the gap A and the 
renormalised Fermi velocity vfr may be extracted; finally Sec. 13.31 presents a fit to a 
similarly-motivated equation of state for A{m,g'^) which along with the assumption of 
hyperscaling permits an estimate of the dynamical critcal exponent z. We summarise 
our findings for the critical parameters in Sec. HI and also attempt to relate our value 

^in [6j the identification of the QCP was restricted to the case {oo, Nfc). 
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for g1 to estimates of ac in the literature. 



2 Formulation and Interpretation of the Model 

Our starting point is a model of relativistic Dirac fermions moving in 2+1 dimensions 
and interacting via an instantaneous Coulomb interaction. In Euclidean metric the 
action is HEIE]: 



dXQd'^x{'lpalQdQ'4)a + VF'4)a^.V'llJa + iVi)al0'^a) + 1^ / c/Xorf'*x((9iF)^ (2) 



a=l 



where e is the electron charge, vp the Fermi velocity, V the electrostatic potential, and 
the 4x4 Dirac matrices satisfy {7^,7;^} = 25^,^, /i = 0,1,2. For monolayer graphene 
the correct number of fermion flavors Nf = 2. The momentum-space propagator for the 
l^-field Di, which couples conserved charge densities ipjo'^ at differing spacetime points, 
is given by 



m , Nj \p\ 



2 



-1 



Di{p)=^ + ^^ (3) 



e2 8 (p2) 



where the first term in brackets on the right hand side is the classical Coulomb interac- 
tion, and the second is the leading quantum correction in the large- A^j limit, describing 
screening due to particle- hole virtual pairs. Note that = Pq + Vp\'p\'^ . The relative 
importance of quantum versus classical effects may be parametrised by the ratio A of 
the two terms in the static limit 0; in SI units 

^ = T7i ^ ^ ^ 

IGeeohvp e 

where e > 1 is the dielectric constant of the underlying substrate. 

For sufficiently large interaction strength the description in terms of massless rel- 
ativistic excitations may be disrupted by condensation of bound fermion-hole pairs in 
the ground state, signalled by an order parameter (V^'?/') ^ 0, with the result that a 
gap appears in the fermion spectrum. Physically this corresponds to a transition from 
a conductor to an insulator; in the language of particle physics the same phenomenon, 
resulting in a dynamical generation of a particle mass, is known as chiral symmetry 
breaking. As this transition occurs at zero temperature, the model predicts a finite 
sequence of quantum critical points (QCPs) whose properties at the critical interaction 
strength Xc{Nf) are sensitive to the value of Nf [6]: the sequence will terminate for Nfc 
(not necessarily integer) defined by X{Nfc) = oo. 

This situation has motivated us to explore a related but distinct model for graphene, 
with action |T0] (in units where vp = ^) 



82=/ I dxod^x 



a=l 



(5) 
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This model resembles the 2+ld Thirring model [HI [12], a four-fermi model known to 
exhibit a sequence of iVj-dependent QCPs as the coupling strength is varied. The 
relation with ([2]) is clarified by inspection of the ^-propagator: 

yg^ 8 (p2)5 J 

Since the quantum correction is identical, models ([2D and ([5]) should yield the same 
physics in the large- A limit, which can be reached as either Nf —* oo or g'^,e'^ oo. 
Indeed, in [lOj we used this property to predict the critical number of fiavors Nfc = 4.8(2) 
above which the model ([2]) remains a semimetal for all g'^. Thus graphene with Nf = 2 
is predicted to be a Mott insulator for sufficiently strong inter-electron coupling. For 
finite Nf the models ([2|) and ([5]) are distinct, although we may still hope they describe 
similar physics for A not too small. 

Let us discuss this point a little further. The principal difference between ([3|) and 
([Hj) occurs at large distances, ie. lim^^oo -Di (r) oc r~^, indicating that the long-range 
Coulomb interaction is not screened, whereas D2 is finite-ranged, being cut off for 
r ^ 0{g'^). It is important to understand whether the modification Di 1— D2 changes the 
physics in any essential way, eg. by defining a model in a different universality class. We 
will be unable to answer this question definitively with the simulation results presented 
here, but note that in the model approach of [5] which predicts dynamical symmetry 
breaking, the relevant momentum range responsible for gap generation is \p\ ^ A/f^, 
implying that it is the short-ranged behaviour of D which governs the properties of the 
QCP. In addition, we note that unlike the "instantaneous" approximation used in that 
work D2 correctly incorporates the po-behaviour of the vacuum polarization function. 

In this paper we will use numerical simulations of a lattice model based opon a 
discretised version of ([Sj) to study the semimetal-insulator transition for the physical 
value Nf = 2. The lattice model is formulated in terms of single-component Grassmann 
fields X; X defined on the sites x of a three-dimensional cubic lattice, by the action 

Siatt = 5^Xx^[(l + Vv^e^''^)Xx+A-(l + '^/.ov^e-^^^-a)x.-A]+m^XxXx. (7) 

X/J, X 

The sign factors rj^^^ = (— l)"^'""* ^-^m-i ensure that in the long-wavelength limit the first 
(antihermitian) term in Siatt describes the Euclidean propagation of Nf = 2 fiavors of 
relativistic fermion described by four- component spinors [14]. The bare fermion mass 
m provides a IR regulator for modes which would otherwise be massless in the limit 
of weak interactions. The hopping terms in Siatt involve the auxiliary boson field Vx 
which is formally defined on the timelike links connecting sites x with x + 0. For further 
details of the relation between the actions ([7[) and ([5]) we refer the reader to [11] and 
[To] H. Because we restrict our consideration to an integer number of fermion fiavors in 

^In particular, for Nf = 2 the action ([7]) can be recast in a "non-compact" form yielding identical 
physics, whose relation to to ([5]) is more manifest. 
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this paper we were able to perform the simulation using a well-established numercial 
method called the hybrid Monte Carlo (HMC) algorithm, which generates equilibrated 
ensembles of field configurations {V} with no systematic bias [T5l [TT]. 

Since our initial paper [IQJ, results from simulation of an alternative lattice approach 
to the graphene model ([2]) have appeared [Zl El [9]. This formulation is based on lat- 
tice gauge theory, in which the degrees of freedom corresponding to the electrostatic 
potential V are formulated on a 3+ld lattice, while the electron degrees of freedom are 
confined to a 2+ld "braneworld" . This in principle gives a more faithful rendering of 
the physics encapsulated in ( P[3i) . The principal result is a prediction for the critical 
coupling corresponding to the semimetal-insulator transition; for Nj = 2 the value 

Af^ = 1.70(2) (8) 

was obtained. This result is intriguing because it lies between the values A ~ 1.25 
expected for graphene on an Si02 substrate, which experiments have shown to be a con- 
ductor, and A ~ 3.4 (Cf. eqn. (jlj)) for freely-suspended graphene, which is accordingly 
predicted to be an insulator. The necessity to store and evolve variables on an extra 
dimension is clearly a computational burden the action ([7j) evades; however for current 
purposes the main advantage we claim for our approach is that it permits a straight- 
forward means to measure the quasiparticle propagator, as presented below in Sec. 13. 2[ 
without the need for gauge fixing. 

Let us finish this section by discussing the behaviour we expect of the model ([7]), and 
how it might relate to physical graphene. In the limit m — the model has a global chiral 
symmetry Xx ^ GxpiaexXx',Xx ^ GxpiaExXx where the sign factor Ex = (—1)^0+^1+^2 
distinguishes odd and even sublattices: the model studied in [3, El [9] has the identical 
symmetry. For Nf flavors the pattern of symmetry breaking expected for the continuum 
models ([21 [5]) is U(2A^/)— i>U(A^/)®U(A^/), whereas away from the continuum limit the 
pattern for ([7]) is U(^)(8>U(^)— >U(^). By analogy with the Thirring model pT] , 
we expect that for large values of the coupling g"^ the symmetry will be dynamically 
broken as signalled by a non- vanishing condensate (xx) = V~^d\n Z/dm 7^ 0, but 
that the symmetry will be restored in a continuous phase transition at some critical 
coupling g'"^. This transition between the two phases defines a UV-stable fixed point 
of the renormalisation group, and the fixed point theory is thus uniquely specified by a 
set of critical exponents - one of the main goals of the paper, presented in Sec. 13.11 is 
to determine both the critical g'"^ and the set of exponents by numerical means. The 
fixed-point theory should describe the low-energy excitations of physical graphene in the 
continuum limit, reached either from the insulating phase as g~^ y g'"^, {xx) ^ 0, or 
from the conducting phase as g~^ \ g~'^, m — >■ 0. As we shall see below, it may be 
possible to relate the value of g~^ to A as defined via (jl]). The applicability of the model 
to graphene, however, rests on the hypothesis that the low-energy excitations and their 
interactions in graphene share its symmetries, and that the physical parameters are such 
that graphene lies within the basin of attraction of the fixed point. Ultimately this must 
be settled by experiment. 
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3 Numerical Results 



Preliminary simulations with Nf = 2 presented in [10] showed evidence for a crossover 
from strong- to weak-coupling behaviour at g'"^ ~ 0.6. Accordingly, we have undertaken 
a refined campaign of simulation on system sizes x Lt with Lg ranging from 16 to 48, 
and Lt ranging from 48 to 84; the bare mass m was varied between 0.0025 and 0.025. 
Because the action ([7]) does not treat space and (Euclidean) time directions equivalently, 
it is useful to explore the consequences of independently varying Lg and Lt: however, 
the most detailed coverage of the axis was obtained on 24^ x 48. 

3.1 Equation of State 

In the vicinity of a second order phase transition the order parameter over a range of 
couplings and small source values m can be described by an equation of state of the 
form 

= At{xxr + B{xx)' + 0{t^xxr'^), (9) 

where g'"^ is the critical couphng, t = g^^ — g^^, is a universal scaling function, S 
the critical exponent describing the order parameter's response at criticality to a small 
applied source m, (5 the exponent governing the scaling of the order parameter for m = 
as t — >■ 0_, and p = 5 — 1/ j3. Order parameter data taken in the thermodynamic limit 
can be fitted to ([9]) to extract g""^, 5 and p. In practice we need to make assumptions 
about the width of the "scaling window" in g~^ and m where the subleading corrections 
in ([9]) can be safely ignored, and we also need to carefully monitor the effects of working 
with finite L^, Lt. 

First let's discuss finite volume effects. Since the model ([7]) has an anisotropic action, 
we cannot a priori exclude the possibilty of correlation lengths in spatial and temporal 
directions diverging with distinct exponents and ut [16]. In previous work we have 
attempted to incoporate this possibility via a correction to the equation of state fit, but 
the complicated nature of the finite volume scaling model made these fits of questionable 
value given the range of simulation volumes available to us. Here we take a more 
pragmatic approach, and compare order parameter data for two different bare masses 
at fixed Lt and varying Lg in Fig. [T^ and vice versa in Fig. [T)d. The plots reveal the 
very different nature of the finite size effects in each case: {xx) rises as Lt is increased, 
corresponding to the zero temperature limit, but falls as the thermodynamic limit Lg 
oo is approached. Moreover, in both cases the effects are greater in the symmetric phase, 
ie. at larger values of g~^. We will proceed by using the observation that in the restricted 
range 0.525 < g''^ < 0.65 the data for m = 0.005 (0.525 < g~^ < 0.70 for m > 0.01) on 
24^ X 48 are free from finite size effects almost within statistical error. 

Table [1] shows sample fits to (Q to 0{t) for order parameter data taken on a 24^ x 48 
lattice, and shows how the fit quality improves as data far from criticality are successively 
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Figure 1: (Color online) {xx) vs. g ^ for m = 0.01 calulated on (a) x 48 lattices; (b) 24^ x Lt 

lattices, showing finite size effects. 



excluded. We also tried excluding low mass points as these are most susceptible to finite 
volume effects. It is comforting that the fitted values of the critical parameters are quite 
stable as the scaling window is so varied. Our preferred fit is the third row of Table [H 
which includes as much data as possible consistent with preserving acceptable fit quality, 
is 

^-2 = 0.609(2); 5 = 2.66(3); p = 1.245(11) {3 = 0.71(2). (10) 

The fitted equation of state is plotted in Fig. [2J 

We also experimented with fits with an extra free parameter modelling an O(t^) 
correction to ([9]); these fits indicated a slightly larger value of g'"^, but despite the 
extra free parameter did not yield appreciably better values. Moreover, the resulting 
equation of state clearly failed to be physically reasonable in the symmetric phase: since 
5 — I < from ( fTOl) . the O(t^) term rapidly becomes numerically dominant here. In 
support of this. Fig. [3] plots order parameter data taken on 24^ x 48 using axes chosen, 
using the critical parameters ffTOj) . to effect data collapse onto the scaling function JF 
which is seen to be linear to very good approximation. 

Finally we examine another probe of the critical point, the ratio of transverse to lon- 
gitudinal susceptibilities Xiixt = d\ia{xx) /dlnm [Hill]. There are two spin-0 particle- 
hole channels with opposite intrinsic parities, which by analogy with mesons in particle 
physics we refer to as a (parity +) and vr (parity — ). In the m ^ limit, in the 
conducting phase the two states are related by a U(l) global chiral symmetry and are 
degenerate; in the insulating phase, by contrast, the symmetry is spontaneously broken 
and the yr-channel therefore contains a massless pole by Goldstone's theorem. Since 
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Figure 2: (Color online) Fit to ([9]) to order parameter data taken on 24^ x 48. The function in the 

m — > limit is also shown. 




- — 1/p 

t<ia> 



Figure 3: (Color online) Plot of m/ {xxY ^s. {g ^ — 5c ^)(xx) using the critical parameters PH)) . 
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fit # 




xVdof 


0.525 < g-'^ < 0.90 (all m) 69 


0.608(2) 2.66(2) 1.252(4) 


6.9 


0.55 < g-'^ < 0.80 (all m) 55 


0.607(2) 2.68(3) 1.261(9) 


4.0 


0.525 < g-^ < 0.65 (m = 0.005) 43 
0.525 < g'^ < 0.675 (m = 0.0075) 
0.525 < g-^ < 0.70 (m > 0.01) 


0.609(2) 2.66(3) 1.245(11) 


2.7 


0.525 < g-^ < 0.70 (m > 0.0075) 37 


0.600(3) 2.80(5) 1.285(14) 


2.3 



Table 1: Various fits to the Equation of State © for data taken on 24^ x 48 



Xe/xt is simply the ratio of the integrated cr-propagator to the integrated vr-propagator, 
we expect it to tend to unity as m — > in the conducting phase, and to zero in the 
insulating phase. Exactly at criticality, however, the ratio is m-independent and takes 
the value 1/S [T7]. Fig. H] plots Xi/Xt vs. m evaluated on a 24^ x 48 lattice, including 
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Figure 4: (Color online) Susceptibility ratio xe/Xt vs. m for various g ^ in the critical region. The 
fitted value of S^^ from pH]) is shown as a horizontal band. 



contributions from diagrams with both connected and disconnected fermion lines as de- 
tailed in prj. The data taken at g^"^ = 0.60,0.625 are approximately m-independent, 
especially for larger m, and bracket the value of obtained from the equation of state 
fit, strengthening our confidence in the values of the critical parameters in (fTOl) . 
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3.2 Quasiparticle Dispersion 

One of the main motivations for the choice of model ([7]) is that since it has no manifest 
gauge symmetry, there is no requirement to fix a gauge in order to define or measure 
a correlation function such as the fermion propagator. This has enabled us to perform 
the first numerical simulation of the quasiparticle excitation spectrum in graphene. 

The fermion excitation spectrum of the model is accessed via analysis of the Euclidean 
timeslice propagator C/(p, t) defined by 

Cfip, t)=Yl (X{0, OMx, t))e-^^^■^ (11) 

X even 

where "even" refers to sites with spatial coordinate x obeying (—1)^^ = (—1)^^ = 1, 
and the components of p take values 27m/ Lg, with n = 0, 1, . . . , Ls/4. This restriction 
improves the signal to noise ratio, and originates in the observation that the action (^^) 
is invariant only under translations by an even number of lattice spacings. The energy 
E{p) is then extracted by a fit of the form 

Cfip,t) = 5(e-^* + e-^(^*-*)), (12) 

where in this case only data with t odd were used, since this yielded the best fits across 
the whole range of g""^ (it can be shown that limm^o C/ = for even t in the conducting 
phase). 

We measured E^p) for p = (pi,0) on 32^ x 48 for g^"^ = 0.55,0.6,0.7,0.8, and 
additionally on 48'^ for g^^ = 0.8, 0.9, using m ranging from 0.005 to 0.03. The resulting 
dispersions for the latter two systems at m = 0.005 are shown in Fig. [51 For small p and 
m the dispersion starts out linear to good approximation, and then flattens out to have 
zero slope at the effective Brillouin zone edge at p = |; this flattening is a discretisation 
artifact with no physical significance. To proceed we parametrise the dispersion relation 
using 

E{p) = Asinh-\^sm^p + M^), (13) 

where for A = 1 and M = m the exact result for non-interacting lattice fermions 
is recovered. Two sample fits are shown in Fig. [51 For small M we can interpret 
£■(0) = A ^ AM as the quasiparticle mass (or gap), and for small p in the limit M — > 
then dE/dp ^ A is the physical Fermi velocity vfr- Results for A and M as functions 
of m are shown in Figs. [6171 

The results for M are broadly consistent with our identification of the critical cou- 
pling. For g^"^ < g~'^ ~ 0.6, Fig. [6l supports limm^oM ^ 0, signalling the generation 
of a gap via spontaneous chiral symmetry breaking. For weaker couplings the data can 
be plausibly extrapolated in the same limit to M = 0, signalling a chirally symmetric, 
conducting phase. Note that throughout the critical region A ^ m, indicating large 
mass renormalisation due to strong interactions even in the symmetric phase. In Sec. 13.31 
below we will present further results for M for a range of m in the critical region. 
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Figure 5: (Color online) Dispersion relation E{p) as measured on a 48^ lattice with m — 0.005. 



Fig. [7] shows that despite some noise in the data the parameter A, and hence the 
physical Fermi velocity vpR, is both m- and ^-independent in the critical region, 
taking a numerical value ~ 0.7. We interpret this as being due to a renormalisation of 
the bare Fermi velocity vp = ^ due to quantum effects. This in principle needs be taken 
into account when we attempt to assign a physical value to the critical coupling g~'^ in 
Sec. m This result is intersting because analytic calculations based on weak-coupling 
and/or large- A^/ predict that vpR > vp 



3.3 Dynamical Critical Exponent 

In this section we take a closer look at correlations in the critical system, both via a high 
statistics study (typically several thousand HMC trajectories) at g'"^ = 0.6, close to the 
critical value reported in (ITU]) , and a refined study of the quasiparticle mass parameter 
M in the critical region. All results are from simulations on 24^ x 48 lattices. Because 
the model ([7]) treats space and time differently, correlation lengths defined in spatial 
and temporal directions can exhibit different critical scaling, leading to two distinct 
exponents defined via \l6l 

^s^\tr-, 6oc|tr^'. (14) 

Our goal in this section is to constrain the value of the dynamical critical exponent 
z = Vt/vg relating spacelike to timelike correlations. 
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Figure 6: (Color online) The fitted parameter M vs. m for various g 
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Fig. [8] shows data for both the order parameter {xx) ^"^^ mass parameter M defined 
by (fT3|) as a function of m on a log-log plot. The linear nature of the plots supports a 
power-law scaling: 

1 

{xx) ^'"^^'i M (X m'^ , (15) 

where 6 and (3 coincide with the definitions implicit in ([9]), and the exponent ut is the 
one relevant for the extraction of spectral properties via (!T2l) from correlations in the 
Euclidean time direction. Least-squares fits (excluding m = 0.025) yield 6 = 2.85(1); 
1^ = 0.38(2). The mismatch between this value for 6 and that of (fTOj) extracted from 
the equation of state is ascribed to the actual value of lying slightly above 0.6, as 
suggested by Fig. HI 

Fig. [9] shows results for the quasiparticle mass parameter M for a range of m and 
g^^ values in the critical region. We have fitted these data with a relation inspired by 
the equation of state ([9]): 

m = AtM-t + BM-t , (16) 

which with M oc .^^^^ understood recovers f|T^ in the limit m — > 0, and is consistent 
with f|T5l) when t = 0. A fit to 33 datapoints with (7^^ fixed by ffTOj) yields 

^ = 2.25(5); ^ = 1.16(6) (17) 

with per degree of freedom of 1.0. 

It is now time to discuss the possible anisotropy at (^"^ = g~'^ in more detail. As 
mentioned above, the ratio z = vt/^s defines the dynamical critical exponent, which is 
an important characteristic of a QCP. In particular, the critical dispersion relation is 
modified to be of the form E oc p^, which has important implications for the stability of 
quasiparticles; energy and momentum conservation make it impossible, in an inelastic 
collision, for a quasiparticle to decay into constituents with smaller E and p if 2; < 1 [B]. 

The results in this section permit an estimate of z via the following indirect argument. 
First, we use the exponent values from f|T7|) and 6,(3 from ffTOj) to estimate 

ut = 0.80(3). (18) 

Next, we use a modified hyperscaling relation [191 [IS! 

i^t + {d~l)us = (19) 

where d = 3 is the number of spacetime dimensions, to estimate 

Us = 0.89(3), (20) 

leading to 

^ = 0.90(5). (21) 

This result is tantalising, since although it hints at 2; < 1 it eliminates neither the value 
z ~ 0.8 based on a leading order large- A^j calculation in the strong-coupling limit [B], 
nor the general result z = 1 claimed for systems at a QCP with d < 4 interacting via a 
Coulomb potential [20] . 
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4 Discussion 



The main result of this paper is that by numerical means we have identified a quantum 
critical point, corresponding to a semimetal-insulator transition, for a model with Nf = 2 
fiavors of Dirac fermion sharing many symmetries with a low energy effective theory of 
monolayer graphene. We have been able to identify critical exponents characterising the 
transition, as summarised in eqn. (fTOj) : the most robust prediction, emerging from a fit 
to the equation of state, and supported by both a calculation of the susceptibility ratio 
Xe/Xt and a direct study of the scaling of the order parameter against m at g'"^ ~ g'"^, 
is that the exponent 6 = 2.66(3). This is significantly different from the value 6 = 5.5(3) 
obtained in the strong-coupling limit at Nf = Nfc = 4.8(2) [10], demonstrating that the 
universality class the model falls into is A^/-dependent. A similar picture has emerged 
from numerical simulations of the 2+ld Thirring model [TTl[T2^j. where it has been shown 
that 6 increases with Nf. Drut and Lahde have reported the same trend from numerical 
simulations of their model with Nf = 0,2,4 [9]. However, their most recent value for 
S{Nf = 2) = 2.26(6) is significantly different from ours, so it remains unclear whether 
the two models lie in the same universality class, or whther the long range interaction 
present in the model of [9J but not here have a decisive effect. 

We have also for the first time presented results for quasiparticle propagation, finding 
evidence for a gap developing spontaneously in the spectrum for g^'^ < g'"^ as m — > 0. 
In addition, analysis of correlations at non-zero momentum has enabled us to roughly 
calculate the renormalisation of the Fermi velocity vp- We reiterate that in our model 
the fermion propagator is uniquely defined and readily calculable; in the original model 
([2]) the presence of a local gauge symmetry makes analysis of quasiparticle propagation 
potentially problematic both theoretically and numerically. 

We have also outlined a method to obtain the dynamical critical exponent z, an 
important characteristic of any QCP, using scaling and hyperscaling arguments. Unfor- 
tunately the inevitable accumulation of errors in such an indirect approach precludes us 
at this stage from conclusively deciding whether z < 1 or not. This issue is of theoret- 
ical interest since there are general arguments to claim z is exactly one for Coulombic 
systems [20] (of course, strictly our model is not in this class). In this respect a more 
direct attempt to extract z via measurements of the quasiparticle dispersion E{p) on 
lattices with a large spatial dimension giving enhanced momentum resolution may prove 
interesting. 

Finally, while our model should be regarded as sharing universal features of graphene 
in the neighbourhood of some putative fixed point, and hence at best able to make 
predictions of critical exponents and dimensionless ratios of low-energy observables, it is 
difficult to resist the temptation to attempt to convert our result g~'^{Nf = 2) = 0.609(2) 
for the critical coupling into a physical prediction. 

First, we must express our result in terms of the continuum model In order 
to do this, we remind the reader that the expression for the propagator is derived 
using a regularisation which respects current conservation; unfortunately the lattice 
regularisation defined by ([7]) is not of this type. The solution, as outlined in [11], is to 
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take the strong-coupling limit of the lattice model not at = but at g~ = g^^, 
which may be identified numerically via the location of a peak in (xx) ^M- The relation 
between the \^-propagator on the lattice and in the continuum is then [ TT] 



Aatt(p; g) = ZD'^{p] gn) = Z 



1 ^ Nf bt 



(22) 



2 

with Z = {1 — |— > 1 and g]^ = Zg^. The extra factor of 2 in the first term in square 

^lim 

brackets results from a careful counting of the staggered fermion degrees of freedom in 

For our graphene model Fig. 1 of Ref. [10] suggests g^^^ = 0.30(2), yielding a renor- 
malised critical coupling g'^j^ ~ 3.26. Now, in order to compare the quantum and classical 
terms in -Diatt to define an effective value of A, a momentum scale p is needed. Since 
the only length scale in the problem is the lattice spacing, a natural (if somewhat ar- 
bitrary) choice is po = 0, IpI = ^; this means that the propagators Di and Z^iatt match 
at a distance of roughly one lattice spacing. The matching condition A = 5'|j;vr/4 yields 
Ac ^ 2.6. 

Since D2 decays faster than Di at large distances, this estimate for Ac is likely to 
be on the high side. We should however note two factors neglected in this simplified 
approach. Firstly, the renormalisation factor Z ^ 2.0 boosts the interaction strength 
of the lattice model; taking proper account of this will have the effect of raising the 
predicted Ac. Secondly, the Fermi velocity vp appearing in ([2]) and implicitly in (jSj) 
is the bare one, whereas presumably it is the renormalised vfr ~ 0.7 vp which has 
the experimental value lO^ms"^. Since A in (jll) is defined in terms of the bare value, 
this correction has the effect of lowering the predicted Ac, although it will also correct 
the A-values calculated for physically-realised cases such as graphene which is either 
freely-suspended or mounted on a substrate of known dielectric constant. In our view 
the uncertainty over the phenomenologically- relevant value of vp must ultimately be 
settled by an ab initio microscopic calculation; moreover, should it prove to be the case 
that z < 1, then the very notion of a universal Fermi velocity becomes ill-defined since 
limp^o ^ diverges. 

Table [2] compares our estimate for the critical interaction strength with both that 
of the alternative simulation of Drut and Lahde [TJ [H [9], with a value predicted using 
a renormalisation group treatment of radiatively-induced four-fermion contact interac- 
tions [2IJ, and with one older and three more recent estimates [H El [221 [231 El] based on 
self-consistent diagrammatic calculations using the model ([2|). Note that is conventional 
in the literature to quote the critical effective fine structure constant Oc = 4Ac/vrA^j; here 
we show both parameters where appropriate. Also listed are estimates for the critical 
number of fiavors Nfc corresponding to the location of the QCP in the strong-coupling 
limit. We argued in [lOj that in this limit our model coincides with (12|) and hence that 
Nfc = 4.8(2) is a robust non-perturbative prediction. 

To give these numbers some meaning recall that A is calculated to be 1.25 for graphene 
on a Si02 substrate, where experimentally it is known to be a conductor, and 3.4 in 
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reference 




Ac 


Nfc 


this work, 


- 


2.6 


4.8(2) 


mm 


1.11(6) 


1.70(8) 


4-6 


m 






2.03 




2.33 


3.66 


2.55 


[22] 


1.13 


1.77 


3.6 


m 


1.16 


1.82 


3.5 


m 


1.62 


2.54 


2.8 



Table 2: Predictions for the critical interaction strength for Nf — 2, and for the critical number of 

flavors Nfc in the strong coupling limit. 

vacuum. Acknowledging the difficulties in obtaining precise numbers reviewed in the 
previous paragraphs, we may nonetheless observe that our simulations lend support to 
the mounting body of theoretical evidence that freely-suspended graphene is an insulator. 
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